Broadening bandgaps in a multi-resonant piezoelectric metamaterial plate via bandgap merging phenomena

Locally resonant metamaterials usually have narrow bandgaps, which significantly limits their applications in realistic engineering environments. In this paper, an optimization method based on the genetic algorithm is proposed to broaden bandgaps in multi-resonant piezoelectric metamaterial through the merging of multiple separated bandgaps. Using the effective medium theory, the equivalent bending stiffness and dispersion relationship of a metamaterial plate are first obtained. Then, the criteria for determining the bandgap ranges for the two cases with and without damping are provided and analyzed. Furthermore, based on the bandgap merging phenomena, an optimization method for widening the bandgap is proposed based on the genetic algorithm. By investigating the bandgap widening effects in cases without and with damping, it is found that, when there is no damping, the bandgap can only be slightly widened; while after introducing damping into the transfer functions, the bandgap can be significantly widened by more than 200%. The bandgap widening effects are verified by comparing with finite element simulation results.

and dispersionless operation.Tsilipakos et al. 9 also proposed a practical, achromatic microwave metasurface for delaying broadband pulses in reflection by fitting five sharply resonant meta-atoms in a subwavelength unit cell, achieving a total bandwidth nearly linearly proportional to the number of resonances.
Piezoelectric materials have also been introduced into the designs of locally resonant metamaterials for their reconfigurability [10][11][12][13] .Airoldi and Ruzzene 10 found in their research that the equivalent stiffness of a piezoelectric beam shunted with resonant circuits exhibits resonance characteristics near the resonance frequency.They proposed that such piezoelectric beams can be regarded as metamaterials, namely piezoelectric metamaterials.Sugino et al. 11 studied one-dimensional locally resonant piezoelectric structures with segmented electrodes under transverse vibration and found that the formation of bandgaps is related to the resonance frequency in the shunting circuit and the electromechanical coupling coefficient of the piezoelectric patches.Yi et al. 12 demonstrated through the analysis of the dynamic properties of piezoelectric metamaterial beams that bandgaps are generated by negative bending stiffness.Similar to passive locally resonant metamaterials, piezoelectric metamaterials also need to address the issue of bandgap widening.For example, Zhu et al. 14 experimentally investigated the dynamic behavior of adaptive piezoelectric metamaterials with negative capacitance and obtained adjustable bandgaps.Sugino et al. 15 proposed a hybrid unit structure made of piezoelectric laminates with segmented electrode pairs and additional mechanical resonators, which can simultaneously exhibit locally resonant bandgaps based on dynamic mass and stiffness.The combination of these two types of bandgaps enhances the overall bandgap.Jian et al. 16 introduced a novel gradient piezoelectric metamaterial beam with parallel resonant circuits, employing a graded strategy for electrode pairs' spatial variation.With proper selection of resistors, the gradient piezoelectric metamaterial beam achieved theoretically the widest attenuation region.Celli et al. 17 proposed that the heterogeneity of resonator properties and the interaction between spatial disorder could lead to the widening of the system's filtering effect, namely widening of the bandgaps, compared to more traditional concepts of spatially ordered rainbow materials.In order to enhance the practical performance of metamaterials, reduce the required number of unit cells, and additional mass, Airoldi and Ruzzene 18 achieved multi-frequency resonance in piezoelectric metamaterials beams and simultaneously generated multiple bandgaps by introducing multifrequency resonant circuits.Furthermore, Yi et al. 19 designed a transfer function to achieve multi-resonance in digital circuits, realizing multi-resonant piezoelectric metamaterials based on self-inductive sensing piezoelectric patches and digital circuits.However, approaches for widening bandgaps in piezoelectric metamaterials are still somewhat limited.
Optimization methods have brought an avenue for the widening of the metamaterial's bandgap.Common optimization schemes include topological optimization and parameter optimization [20][21][22][23][24][25] .Regarding passive elastic metamaterials, there is a considerable amount of research, for instance, Wu et al. 22 proposed a genetic algorithm approach to design non-periodic unit structures, achieving broadband wave attenuation.Meng et al. 23 established a theoretical model for connecting multiple local resonance bandgaps in two-phase composite materials with damping and obtained ultra-wide bandgaps through the proposed optimization procedure.Chen et al. 24 employed an enhanced genetic algorithm to design square lattice structures, achieving wide and multiple bandgaps in the low-frequency range through optimized filler material distribution.For piezoelectric metamaterial, Jian et al. 25 introduced a non-uniform piezoelectric metamaterial beam, where the shunting circuit parameters were optimized using an adaptive genetic algorithm to tailor the vibration attenuation region.Gao and Wang 26 proposed a genetic algorithm to optimize the circuit parameters of LR parallel circuits in each unit of the metamaterial beam, which can couple multiple local resonance bandgaps to Bragg bandgaps to widen the frequency range of vibration suppression.In general, it remains an unresolved issue to widen bandgaps and achieve broadband vibration suppression in piezoelectric metamaterial plate through optimization methods.
This paper proposes a method to widen the bandgaps of multi-resonant piezoelectric metamaterials based on frequency merging phenomena.Firstly, a theoretical model of a piezo-metamaterial plate is obtained based on the effective medium theory, and criteria for determining the bandgap range in the metamaterial are provided by considering scenarios with and without damping.Subsequently, the bandgap merging phenomena between separated multi-bandgaps are analyzed.Then, based on the genetic algorithm and considering the bandgap merging phenomena, an optimization design method for widening the bandgaps is proposed.The bandgap widening effects by using the optimization method are studied separately for transfer functions with and without damping.Finally, numerical simulations are conducted to verify the effectiveness of widening bandgaps in the metamaterials.

Multi-resonant piezoelectric metamaterial
Metamaterial plate physical model Figure 1 depicts a schematic diagram of the piezoelectric metamaterials plate and its unit.Each unit consists of a base plate, two piezoelectric patches attached to the upper and lower surfaces of the plate structure respectively, and digital synthesized impedance circuit connected to them.The metamaterial plate is composed of units arranged periodically in the plane.The polarization direction of the piezoelectric patches is aligned with the z-axis, and the surface electrode in contact with the plate is grounded.The digital synthesized impedance circuit includes voltage scaling circuit, input voltage biasing circuit, controller, output voltage biasing circuit, and voltage-controlled current source 27 , the overall effect of the circuit can be represented by the transfer function G.The material parameters and geometric dimensions of the base plate and piezoelectric patches in a metamaterial unit are shown in Table 1.
www.nature.com/scientificreports/Bandgap determination criterion Firstly, the equivalent parameters of the piezoelectric patch with shunting circuits are studied.The piezoelectric patch can be considered as isotropic thin plate.Based on the plane stress assumption and the sub-wavelength assumption, their equivalent in-plane Young's modulus and in-plane Poisson's ratio are 28 : where E sc p and ν sc p are the in-plane Young's modulus and Poisson's ratio of the piezoelectric patch in the shortcircuit state respectively, k 31 = d 31 E sc p /ε σ 3 is the piezoelectric coupling coefficient, C T p = A p ε σ 3 /h p is the intrin- sic capacitance of the piezoelectric patch under steady stress, C s p = C T p 1 − k 2 31 is the intrinsic capacitance of the piezoelectric patch under steady strain, s is the Laplace parameter, d 31 is the piezoelectric coupling constant of the patch, ε σ 3 is the dielectric constant under steady stress, A p is the area of the patch, and h p is the thickness.In the unit, according to classical laminated plate theory 29 , the equivalent bending stiffness of the laminated plate covered by the piezoelectric patch can be expressed as: (1)  where is the bending stiffness of the elastic plate, E b and ν b are the Young's modulus and Poisson's ratio of the elastic plate, respectively.
Due to the different dimensions of the piezoelectric patch and the elastic plate, the equivalent parameters of the entire unit are calculated using the effective medium method 28 : w h e re χ = l 2 p /l 2 b i s t h e c ove r a ge r at i o of t h e pi e z o e l e c t r i c p atch i n t h e u n it , is the equivalent bending stiffness of the laminated plate covered by the piezoelectric patch when the piezoelectric patch is short-circuited, ρ b and ρ p are respectively the densities of the elastic plate and the piezoelectric patch.The transfer function G is the multi-resonant transfer function designed by the authors in 19 .It makes the equivalent Young's modulus E p of the piezoelectric patch resonate at multiple frequencies, and multiple bandgaps are generated in the metamaterial.The expression of the transfer function is: β i , ω z,i , and ω p,i represents the damping, zeros, and poles in the transfer function G respectively, and they correspond one-to-one, with their quantity determined by the number of poles n.The poles can be expressed as ω p,i = 2πf i , f i corresponding to the resonance frequencies of the poles.Then, n bandgaps at the n resonance frequencies corresponding to the poles are produced at the metamaterial.The transfer function should satisfy the Nyquist stability criterion, i.e., its poles should all be located in the left half-plane of the complex plane, and the relationship between poles and zeros is obtained 19 : The frequency dispersion relationship of the metamaterial can be expressed by the equivalent properties of it.Here, we consider the dispersion relation in the x-direction of the metamaterial.The dispersion relation of the metamaterial with two poles configured in the transfer function is shown in Fig. 2, where Fig. 2a and b respectively depict the variation curves of the real and imaginary parts of the wavenumber with frequency.The dispersion relation is described using complex wavenumbers: www.nature.com/scientificreports/Assume a harmonic wave solution u = Ae j(kx−ωt) , from Eq. ( 8), we can obtain: This indicates, the real part of the wavenumber represents the phase change of the wave, while the imaginary part of the wavenumber represents the attenuation of the amplitude in space.
Combining Eqs. ( 4), (5), and ( 8), it can be observed that when the material and geometric dimensions of the metamaterial are determined, the wavenumber value is only related to the transfer function G.According to the expression of the transfer function ( 6), the parameters affecting wavenumber are n , β i , ω z,i , ω p,i .Pole ω p,i and zero ω z,i affect the location where the poles are generated, while damping β i directly affects the wave propagation effect.
In Fig. 2, there are three curves each for the real and imaginary parts of the wavenumber dispersion curves, with different damping configurations.In Fig. 2a, as the damping increases, the resonance effect of the real part of the wavenumber gradually weakens.In Fig. 2b, when the damping is zero, there is a clear range of zero values and non-zero values for the imaginary part of the wavenumber, where the non-zero range represents the bandgap range.And if the damping is non-zero, there is no longer a strict zero-value range for the imaginary part of the wavenumber.Therefore, in these cases, we define the regions where the imaginary part of the wavenumber Im(k) is larger than a certain value k th as bandgaps.We require that the vibration decreases by at least − 10 dB after passing through four unit cells.According to Eq. ( 9), 20 log x 1) = −10 (subscripts 1 and 2 represent two points in the direction of vibration transmission, and the distance between them is four unit lengths), which gives the imaginary part of the wavenumber to be greater than 8.5534.Therefore, we have chosen 10 as the bandgap criterion value here, i.e., k th = 10.Results in Fig. 2 (and other references 19 ) also demonstrate that, in terms of wave attenuation, damping has positive effect between the two poles but has negative effect at and near the vicinity of the two poles.Therefore, increasing damping blindly is not advisable.Instead, damping should be optimized to obtain the best attenuation characteristics.

Bandgap merging phenomenon
When multiple poles are configured in the multi-resonant transfer function, piezoelectric metamaterials can generate multiple bandgaps.It has been found in studies that there is a coupling interaction between these multiple bandgaps.When the bandgaps are close enough, they merge into one, and the merged bandgap is wider than the original multiple individual bandgaps 19 .
Taking the case of two poles configured in the transfer function as an example, the phenomenon of bandgap merging is illustrated in Fig. 3.When the damping in the transfer function is zero, the bandgap merging is dominated by the distance between poles, as shown in Fig. 3a.If the distance between poles is less than or equal to α * = 0.028, two independent bandgaps merge, and the merged one is wider than the original two separate bandgaps.When there is damping in the transfer function, besides the influence of the distance between poles, the range of the bandgap is mainly affected by the damping magnitude, as shown in Fig. 3b.The distance between poles α in the transfer function is set to 0.04, and the equal damping corresponding to the two poles is taken as the independent variable, namely β 1 = β 2 = β .It can be observed that, without changing the distance between poles, as the damping increases, the widths of the two independent bandgaps increase, and bandgap merging occurs when the damping equaling a certain value.It should be noted that after bandgap merging, increasing the damping continuously leads to a decrease in the bandgap width.From Fig. 3, it can be concluded that widening the bandgap can be achieved by merging multiple bandgaps in multi-resonant piezoelectric metamaterials.www.nature.com/scientificreports/When the number of poles in the transfer function exceeds 2, the situation of bandgap merging becomes more complex.This paper considers using optimization algorithms to find the optimal effect of bandgap merging when configuring multiple poles in the transfer function.

Optimization method
Due to the different methods of determining the bandgap range, it is necessary to discuss separately for β i = 0 and β i > 0 in the optimization calculation, with their corresponding bandgap criteria being Im(k) > 0 and Im(k) > k th , respectively.
The optimization aims to merge multiple independent bandgaps generated by multiple poles into a wider bandgap.Since the merged bandgap is wider than the individual bandgaps before merging, the optimization objective is to maximize the length of the frequency range of the widest bandgap among all the bandgaps in the multi-pole metamaterial.After determining the parameters of the transfer function, the imaginary part of the wavenumber values can be obtained from Formula (8).From these values, the frequency range of the bandgap is identified, and the left boundary frequency of the widest bandgap is denoted as ω 1 , while ω 2 is for the right boundary.Therefore, the objective function for optimization is formulated as follows: where the variable parameters are n , β i , ω z,i , and ω p,i .To improve optimization efficiency, it is necessary to con- sider additional constraints to limit the design space.
In the case of two-pole transfer function, the bandgap merging in metamaterials occurs only when the separation between poles is less than a critical value α * .However, in situations involving multiple poles (n > 2), due to coupling effects 19 , the critical value for pairwise merging of bandgaps may be smaller than α * .Therefore, choosing α * to limit the range of pole value search is rational.
For the cases with damping, the bandgap merging phenomenon depends on the value of β and the distance between two poles.But, if the distance between two poles is too large, no matter what value the β is chosen, the two bandgaps cannot be merged.Therefore, in our optimization procedure, to make the algorithm more efficient, we restrict the distance of two poles to be less than a certain value.In our studies, we choose 2α * , and the results demonstrate that we can find all the optimized parameters using this value.In other studies with different material and geometry parameters, one may need to increase this value if the optimization algorithm fails.A metric coefficient p is introduced for α * , namely p = 1 when β i = 0 , and p = 2 when β i > 0 .When considering damp- ing in the transfer function, an appropriate damping factor β should be selected to adjust the bandgap effect.However, the damping coefficient should not be too large to avoid weakening the resonance effect of the shunt circuit.Here, the search range for β is chosen to be less than 0.1.The searching range of β must be enlarged if the algorithm fails to find the optimized parameters.The constraints on variables ω p,i and β i limit the design space to: where ω 0 is the target angular frequency, ω 0 = 2πf 0 , the metamaterial generates a bandgap at frequency f 0 to suppress structural vibrations.During the optimization process, ω p,1 is fixed as the target frequency.For conveni- ence in the optimization process, the pole values are normalized, by setting ω ′ p,i = ω p,i /ω p,1 .In this paper, genetic algorithms are utilized to solve optimization problems with a specified objective J.The principle of genetic algorithms mimics the natural evolutionary process, where individuals better adapted to the environment have higher chances of producing offspring for the next generation, gradually evolving towards superior solutions.The evolution process begins with a randomly generated initial population containing Q individuals, with each individual represented by a vector ω p,2 , . . ., ω p,n , ω z,1 , . . ., ω z,n , β 1 , . . ., β i (which becomes ω p,2 , . . ., ω p,n , ω z,1 , . . ., ω z,n when damping terms are disregarded), where the parameters in the vector are subject to constraints (7) and (10).The optimization objective function ( 9) is directly employed as the fitness function to evaluate the quality of individuals.Based on the obtained fitness values, three operators-selection, crossover, and mutation-are applied to individuals in the population to generate better offspring.Subsequently, a new generation is produced, and the aforementioned process is repeated.Through successive iterations, when further improvements become unattainable, the optimal solution is obtained.Figure 4 illustrates the flowchart of the optimization algorithm design.

The damping term in the transfer function is zero
The optimization strategy above is used to optimize the case of β i = 0 .In this case, the number of variables is (2n − 1) .Taking the case with three poles, the constraint on variable ω p,i is set to 1 < ω p,i ≤ [1 + (i − 1) • 0.028] (i = 2, 3) .The population size is set to Q = 50.Figure 5a illustrates the evolution of fitness values during the optimization process.Here are the optimal parameter values obtained through optimization when the number of poles n is 3: To clearly illustrate the evolution of bandgap merging, the dispersion curves of metamaterial unit cells are used to represent the process, with parameters chosen from the initial generation, intermediate generations, and optimal results, as shown in Fig. 5b-d.It can be observed that in the initial generation, Fig. 5b, there are three independent bandgaps, while in the optimal result, Fig. 5d, only one bandgap remains, formed by the merging of the three independent bandgaps, and its width is greater than that of the individual bandgaps.
As is well known, genetic algorithm is time-consuming, with much of the computational resources spent on repeatedly calculating fitness values.Therefore, the convergence time of the algorithm strongly depends on the population size ( (2n − 1) × Q ) and the complexity of the fitness function.The larger the number of poles n, the more time is required for computation.Therefore, here we first consider cases with 10 poles for calculation.The optimization results are shown in Table 2.It can be seen that the width of the merged bandgap after combining multiple poles is slightly wider than that of a single bandgap.
The metamaterial generates bandgaps in the structure through locally resonant units, which achieve structural vibration control by blocking the propagation of elastic waves within the corresponding frequency range.Then, the elastic wave transmission isolation characteristic of the bandgap is applied, and the optimized metamaterial parameters are selected to simulate the vibration transmission isolation effect in the metamaterial plate in the frequency domain, so as to verify the optimization results.
A simulation model of a piezoelectric metamaterials plate is established in COMSOL6.0, as shown in Fig. 5e.The plate is composed of a 12 × 12 array of metamaterial units illustrated in Fig. 1b, with material parameters and geometric dimensions listed in Table 1.In the metamaterial plate, piezoelectric patches are attached to the upper and lower surfaces of the substrate plate.Add Solid Mechanics and Electrostatics as two physical fields, and couple them to form the piezoelectric effect through multi-physics coupling.Apply free boundaries to the metamaterial plate in the solid mechanics module.In the charge module, ground the surface of the piezoelectric element adhered to the substrate plate.Adding a weak boundary contribution on the outer surface of the piezoelectric patch to connect the circuit with compliance G to the piezoelectric patch in the COMSOL model.Weak contributions are functionalities used internally in COMSOL for applying built-in domains and boundary conditions.Here, the specific weak expression is given by test(V )•V iA p ω • G .It represents the boundary condition of adding charges on the outer surface of the piezoelectric patch.Here, A p is the area of the surface of a piezo-patch, G is the transfer function in Eq. ( 6), V is the voltage of the DOF, test(•) is an embedded function in COMSOL.
The central position of the simulated plate model is denoted as point A, which a unit harmonic point load F is applied to simulate the disturbance applied to the structure, within the solid mechanics module.Point B is Taking the example of three-pole, the theoretical calculation of the normalized bandgap range is [0.9971, 1.0591] .Figure 5f shows the frequency response curve of point B under the effect of the transfer function with three-pole (solid line), while the shaded region represents the theoretical bandgap range, and the dashed lines denote the normalized frequency values of the simulated bandgap boundaries.The displacement response expresses the vertical displacement |w| of the plate caused by its bending stiffness.Three frequencies are selected near the theoretical bandgap boundaries to observe the corresponding vibrational modes, as shown in Fig. 5g.From Fig. 5f, it can be observed that the frequency response function (FRF) of the numerical results exhibits a trough within the theoretical bandgap range, indicating a significant reduction in displacement at the measurement point.It indicates that this interval represents the bandgap range, as the bandgap blocks the propagation of elastic waves within the plate, preventing the vibration generated by the central excitation of the plate from propagating to the measurement point in the far field.The vibrational mode plots in Fig. 5g provide a more intuitive demonstration of the bandgap effect.At 0.997, the metamaterial plate undergoes significant deformation, while at 0.9972, the deformation is localized at the center of the plate, indicating that the frequency is within the bandgap, and the bandgap blocks the elastic wave propagation generated by excitation.The theoretical and numerical bandgap boundary frequencies are in good agreement, with the left boundary frequency matching well.At the right boundary, beyond the theoretical bandgap boundary frequency of 1.059, the deformation of the metamaterial plate increases, suggesting that 1.059 can be considered as the right boundary frequency of the bandgap.Therefore, the theoretically calculated bandgap range is in good agreement with the numerical simulation results.
In addition, the results of n = 6, 9 are presented in Fig. 6, with three plots corresponding to each case.These plots depict the evolution of fitness values during the optimization process, the dispersion curves of the metamaterial under the optimal parameters, and the frequency response curves.In the frequency response curves (c) and (f), the theoretical bandgap range is also represented by shaded regions, while the simulated bandgap boundaries are depicted by dashed lines with normalized frequency values.The bandgap determination method is the same as that of the three-pole transfer function.For n = 6 , the theoretical bandgap range is [0.9975, 1.061] , while the simulated bandgap is [0.9975, 1.061] .For n = 9 , the theoretical bandgap range is [1.0, 1.0624] , and the simulated bandgap is [1.0, 1.062] .The results show that the theoretical and simulated bandgap ranges are in good agreement.
The above analysis results show that the merging of multiple bandgaps can indeed widen the bandgap, but the effect is not very significant.Then, the reason why this method cannot significantly widen the bandgap is analyzed in detail.
The inequality (7) provides a constraint relationship between zeros and poles, which can be reformulated as follows: The range of value for the ratio in Eq. ( 13) narrows as the number of poles n increases.When considering no damping, the single-pole transfer function has only two variables: pole ω p and zero ω z .Analyzing the variation of bandgap width with the zero-pole ratio ω z /ω p , as shown in Fig. 7, it can be observed that the bandgap width decreases monotonically with the increase in the zero-pole ratio.Equation (13) indicates that as the number of poles n increases, the zero-pole ratio increases compared to the left boundary of the value range.This implies that the range of the zero-pole ratio becomes smaller and smaller.As shown in Fig. 7, the corresponding bandgap width also becomes narrower.The width of each bandgap in metamaterials decreases as the number of poles n increases.This also leads to the fact that increasing n does not significantly increase the width of the merged bandgap.

Taking into account the damping effect in the transfer function
Using the optimization strategy mentioned above, optimization calculations are performed for the case of β i > 0 .The number of variables is (3n − 1) .Taking n = 3 as an example, the constraints for variable ω p,i are set as 1 < ω p,i ≤ [1 + (i − 1) • 0.048](i = 2, 3) , and the constraints for variable β i are set as 0 < β i < 0.1 .The number of individuals is set to Q = 100.Below are the optimal parameter values obtained from the optimization when n = 3:   Figure 8a depicts the evolution of fitness values during optimization and the corresponding dispersion curves for the final optimized parameters.In the subplot, the solid line represents the dispersion curve, while the dashed line indicates k th = 10.
Here, the case with ten poles is also considered for calculation.The optimization results are shown in Table 3: Comparing optimization results with and without damping.It can be observed that when considering the damping effect, the width of the merged bandgap of the multi-pole metamaterial increases with the number of poles n, but the rate of increase gradually slows down as n increases.Under the same n, the presence of damping significantly widens the bandgap.For instance, in the case of a damped ten-pole metamaterial, the bandgap is 200% wider than the without damping bandgap.
Similarly, the above optimization results are verified through simulations in COMSOL6.0.Due to the presence of damping, it is not possible to directly determine the bandgap position from the mode diagram of the metamaterial plate.Therefore, we constrain the frequency band where the imaginary part of the wavenumber is greater than 10 to indicate the bandgap.Here, the bandgap position is verified through the wavenumber values obtained from simulation results.Using the same approach as the previous section, the simulation model is established.However, in this case, wavenumber calculation is required to consider the propagation of elastic waves.Therefore, perfectly matched layers are added around the metamaterial plate to absorb reflected waves, as indicated by the dark region surrounding the plate in Fig. 8b.Specifically, establishing a fictitious domain in the Definition, adding a perfect matching layer, and selecting the domain as all substrates and piezoelectric patches in the dark-colored regions.The metamaterial plate, excluding the perfectly matched layers, consists of a 10 × 10 array of metamaterial units, with excitation applied at point A and measurement taken at point B.
In the simulation model, a unit harmonic load is applied at point A in the center of the metamaterial plate, generating elastic waves in it.From the wave solution in Eq. ( 9), the amplitude ratio between point B and point A can be obtained as: , where l b represents the geometric dimensions of the base structure in the metamaterial unit.This leads to the requirements for the amplitude ratio between point B and point A corresponding to the bandgap range: In the COMSOL Definition module, set up two point integrations (intop), selecting points A and B. Define the variable V 1 = |u B |/|u A | to compute the displacement ratio between the two points.After conducting frequency domain analysis, we can plot the curve in the results.
Figure 8c investigates the metamaterial with a transfer function featuring three-pole.The solid line represents the variation of the displacement amplitude ratio between point B and point A with frequency in the simulation results, while the dashed line indicates the required ratio.It can be observed that the range satisfying the above ratio requirement is [0.993, 1.083] , and the theoretically obtained bandgap is [0.9934, 1.0881] (Indicated in the figure by shaded regions), the left boundaries of both are essentially the same, while the simulated results on the right boundary are slightly smaller.However, the error between the simulated bandgap width and the theoretical bandgap width is less than 5%, so the simulated results can be considered reasonable.
Additionally, the results for n = 6 and n = 9 are presented in Fig. 9, including the evolution of fitness values during optimization, the dispersion curves of the metamaterial obtained under the optimal parameters, and the variation of the amplitude ratio between two points in the corresponding simulation results.Dashed lines are used in the amplitude ratio curves (c) and (f) to indicate the bandgap ranges, applying the method described by Eq. ( 16) for determination.The theoretical bandgap range for n = 6 is [0.9971, 1.1103] , while the simulated bandgap is [0.998, 1.105] .Similarly, for n = 9, the theoretical bandgap range is [0.9982, 1.1175] , and the simulated bandgap is [0.999, 1.113] .The left bandgaps in both cases are essentially the same, while the simulated results for the right bandgap are slightly smaller but within the acceptable error range.Therefore, the simulated results can be considered reasonable.
In the previous analysis, the frequency range where the imaginary part of the wave vector is greater than 10 is defined as the bandgap.Here, we discuss the influence of wavenumber selection on the bandgap range for n = 3 by considering other wavenumber values.Optimization calculations are performed for k th = 5 and k th = 15, and the results are shown in Fig. 10.Combining with Fig. 8, it can be observed that when a smaller k th is chosen for optimization, a wider bandgap range is obtained; whereas, with a larger k th , the resulting bandgap range is narrower.In real applications, k th is determined according to the desired wave attenuation level.For example, if we want at least 50% reduce of the amplitude of a wave when it travels from point A to B, we can find the value of k th by k th ≥ − ln |u B | |u A | /(l AB ) = − ln (0.5)/(l AB ) , l AB is the distance between point A and B. The optimization results obtained under three different wavenumber constraints are compared through simulations in COMSOL.Using the infinite model depicted in Fig. 8b, a unit point load is applied at the center of the metamaterial plate, and the frequency responses measured at the same position on the periphery are shown in Fig. 10e.The imaginary part of the wavenumber represents the decay of the amplitude in space.From the figure, it can be observed that within the narrowest band corresponding to k th = 15 (i.e., the range [0.9964, 1.0615]), the response curve represented by the bottom blue line, which corresponds to the constraint of k th = 15, exhibits the smallest response amplitude and the fastest decay.Conversely, the response curve for k th = 5, indicating the smallest wavenumber constraint, exhibits the largest response amplitude and the slowest decay.However, from a wider frequency range, the response curve for k th = 5 shows attenuation over a larger range, while the response www.nature.com/scientificreports/curve for k th = 15 has the smallest impact range.Reducing the constraint on the wavenumber can increase the range of attenuation, but it also weakens the attenuation effect.

Conclusion
This paper presents a method for widening the bandgap of multi-resonant piezoelectric metamaterials.The equivalent properties of the metamaterial plate are derived, and the dispersion relation represented by complex wavenumbers is obtained.By analyzing the influence of damping effects on the dispersion relation in the transfer function, criteria for determining the bandgap range with and without damping in the transfer function are provided.By discussing the phenomenon of bandgap merging in metamaterials, an optimization design method for widening the bandgap in metamaterials based on genetic algorithms is proposed.Optimization designs are conducted for both with and without damping multi-pole metamaterials, and in both cases, the optimized bandgaps are wider than the previous independent ones.Without damping, the bandgap can be slightly widened, while with damping, the bandgap frequency range can be significantly expanded by over 200% through bandgap merging.This validates the effectiveness of the method for widening the bandgap in piezoelectric metamaterials.The proposed method for widening the bandgap in metamaterials can achieve broad vibration suppression capabilities at the target frequency and can actively tune the frequency range of vibration suppression.

Figure 1 .
Figure 1.Sketch of the designed piezoelectric metamaterial plate: (a) metamaterial plate; (b) top view of unit; (c) side view of unit; (d) external digital circuit of unit 27 .

Figure 2 .
Figure 2. Dispersion curves corresponding to a two-pole transfer function with different level of damping.

Figure 4 .Figure 5 .
Figure 4. Flow chart of optimization algorithm to calculate optimal bandgap.

Figure 7 .
Figure 7. Bandgap width changes with the zero-pole ratio.

Figure 8 .
Figure 8. Bandgap merging results of damped three-pole metamaterial: (a) theoretical calculation results: the main plot shows the convergence of the optimized fitness value evolution, and the inset displays the dispersion curve under optimized results with parameters given in Eq. (14); (b) the metamaterial plate model used for simulation with damping; (c) the ratio of points B to A obtained from the simulation as a function of frequency.

Figure 10 .
Figure 10.Detailed explanation of the of k th , using the three-pole metamaterial as an example: (a, b) show the convergence plot of the optimized fitness value and the dispersion curve of the metamaterial for k th = 5, respectively; (c, d) show the results for k th = 15; (e) summarizes the displacement responses under optimized results for k th = 5, 10, and 15.

Table 1 .
Geometric and material parameters.